%{
                                                        histProcess.m V1
                                                     Written by Rajeev Kant
                                                        Sept. 4th, 2020
                                    Coulombe Lab for Cardiovascular Regenerative Engineering
                            Center for Biomedical Engineering, Brown Univeristy, Providence, RI, 02903

This function calculates the orientation angle of nuclei relative to their
location of outgrowth from the aortic ring.

Inputs:
1. region: IF or R region for determining original ring location
2. stats: regionprops stats of subROI to be processed with majoraxislength,
minoraxislength, and orientation fields
3. ROIpos: Position information of subROI being run determined in getROIs.m
4. radRing: radius of ring determined in ARA_Analyze.m
5. min: minimum aspect ratio for a particle to be considered aligned
6. max: max aspect ratio to filter out overly aligned particles (debris,
under/oversegmented particles)
7. segIm: subROI to be analyzed
6. regName: Character array for region ('IF' or 'RM') for naming
9. num: subROI for naming
10. path: file path

Outputs:
1. stats: Saves histProcess outputs into the input stats file
consisting of:
 1a. AR: calculated aspect ratio of particle
 1b. orHist: Relative orientation of particle


orientation fields) and:
1. Add an aspect ratio field to the regionprops struct
2. Filter out particles with overly small or large AR
3. orHist: orientations for all other particles
4. histSum: Summary metrics with mean, STD, and SEM of orHist
length(orHist) / length(stats) will provide % of polarized cells.


Original publication for citation: 
Kant, R. J., Bare, C. F., and Coulombe, K. L. K. "Tissues with 
patterned vessels or protein release induce vascular chemotaxis
in an in vitro platform", Tissue Eng. Part A, 2021.
doi:10.1089/ten.TEA.2020.0269

%}
function [stats] = histProcess (region, stats, ROIpos, radRing, min, max, segIm, regName, num, path)
%% Calculate aspect ratio and add to stats struct

AR = ([stats.MajorAxisLength]./[stats.MinorAxisLength]).';
AR = num2cell(AR);
[stats(:).AR] = AR{:};

% Define min and max values for exclusion: default 1.1 and 5
ARmin = min;
ARmax = max;

%% Calculate transformation angle based on location of subROI

if strcmp(regName, 'IF')
    x1 = size(region,2) + (radRing*7/8); % Initial coordinates of center of image, accounting for ring diameter and IF or RM region
else
    x1 = -(radRing*7/8);
end
y1 = size(region,1) / 2;
x2 = ROIpos(1,1) + (ROIpos(1,3) / 2); % Coordinates for middle of subROI
y2 = ROIpos(1,2) + (ROIpos(1,4) / 2);
X = [x1 x2];
Y = [y1 y2];
theta = atand((y2-y1)/(x2-x1)); % Angle between ring origin and ROI

%% Calculate transformed angles and add to stats

orHist = NaN(size(stats));

% If aspect ratio of particle i is between thresholds, add to orHist and
% transform, otherwise set to NaN
for i = 1:length(stats)
    if stats(i).AR > ARmin && stats(i).AR < ARmax
        orHist(i,1) = stats(i).Orientation + theta; % Transform orHist orientation
        
        % Restrict Nuclei Angle to +/- 90 degrees (since we don't know directionality of nuclei orientation)
        if orHist(i) > 90
            orHist(i) = orHist(i) - 180;
        elseif orHist(i) < -90
            orHist(i) = orHist(i) + 180;
        end
    end
end

orHist=num2cell(orHist);
[stats(:).orHist] = orHist{:}; % Add to stats

%% Draw plot with arrows on aligned nuclei

s = size(stats, 1);
xcoords = zeros(s, 2);
ycoords = zeros(s, 2);
vectx = zeros(s, 1);
vecty = zeros(s, 1);


for i = 1:s
    if isnan(stats(i).orHist) % Filter out Nuclei that are not aliged (AR beyond set boundaries)
    else % Calculate vectors to overlay on each particle
        hlen = stats(i).MajorAxisLength/2;
        xCentre = stats(i).Centroid(1);
        yCentre = stats(i).Centroid(2);
        cosOrient = cosd(stats(i).Orientation);
        sinOrient = sind(stats(i).Orientation);
        xcoords(i, :) = xCentre + hlen * [cosOrient -cosOrient];
        ycoords(i, :) = yCentre + hlen * [-sinOrient sinOrient];
        vectx(i,1) = xcoords(i,2) - xcoords(i,1);
        vecty(i,1) = ycoords(i,2) - ycoords(i,1);
    end
end

imshow(segIm);
hold on;
quiver(xcoords(:,1), ycoords(:,1), vectx, vecty, 0, 'Color', 'r')
iptsetpref('ImshowBorder','tight');
saveas(gcf,fullfile(path, strcat(regName,'_hist',num,'.tif')),'tiff');
hold off;
close all;  

end % End of histProcess.m